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Relativistic MHD simulations of pulsar bow-shock nebulae 
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Abstract. Pulsar bow-shock nebulae are a class of pulsar wind nebulae (PWNe) that form when the pulsar wind is confined 
by the ram pressure of the ambient medium, and are usually associated with old pulsars, that have already emerged from the 
progenitor Supernova Remnant (SNR). Until a few years ago these nebulae were mainly observed as Ha sources; recently, also 
non-thermal emission has been detected. This is the signature of accelerated particles gyrating in a magnetic field. In the same 
way as Ka radiation is a tool for studying the layer of shocked Interstellar Medium (ISM), the non-thermal radiation might 
be used to infer the properties of the shocked pulsar wind. However theoretical and numerical models have been presented 
so far only in the hydrodynamical (HD) regime, while in order to properly model the internal flow structure and the emission 
properties of these nebulae a magnetohydrodynamical (MHD) treatment is required. We present here relativistic MHD (RMHD) 
axisymmetric simulations of pulsar wind bow-shock nebulae. The structure and fluid dynamics of such objects is investigated 
for various values of the pulsar wind magnetization. Simulated synchrotron maps are computed and comparison of the emission 
pattern with observations is discussed. 
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1. Introduction 

Since the publication of the earliest sample of pulsars (Gunn 
& Ostrikerri970), it has become apparent that many of them 
are characterized by a high spatial velocity. With the increase 
of the number of known pulsars, and the improvement in the 
determination of their distance and position, also the largest 
inferred velocity for this class of objects has increased and is 
presently of order 1000 km s'^Cordes & Chernoff 1998 ). The 
velocity dispersion is also noticeable, since it seems to exceed 
that of any other galactic stellar population. An analysis by 
Hansen & Phinney ( 1 997J estimates a dispersion of order 200 
kms"\ while more recent studies (Arzoumanian et al. 2002) 
show strong evidence for a double peaked distribution, with 
characteristic peak velocities of ^ 90 km s"^and ^ 500 km s"^ . 
In spite of the large uncertainties that still affect it, the latter 
result leads to estimate that up to 10% of the pulsars younger 
than 20 kyr are located outside the SNR left by their progenitor. 

Pulsars are known to be source of ultrarelativistic mag- 
netized winds (see e.g. Michel & Li 1999), mainly made of 
pairs and a toroidal magnetic field. Given the high Lorentz fac- 
tor with which the wind is expected to expand, in the range 
10"^- 10^ (Kennel & Coroniti 1984 ), its interaction with the am- 
bient medium creates a region of hot shocked plasma, that may 
be revealed as a source of synchrotron and Inverse Compton ra- 
diation. Emission of this nature is observed in plerions, where 
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a young pulsar is still inside the shell of the progenitor SNR: 
the best known examples are the Crab and Vela nebulae. 

In the case of a pulsar moving through the ISM, the in- 
teraction of its wind with the ambient medium gives rise to 
a cometary like nebula (Bucciantini & Bandiera 120011) . simi- 
lar to what is expected for the solar wind heliosphere (Zank 
11999ft . Typical velocities of pulsars are much higher than the 
sound speed in the ISM, which ranges from ~ 10 kms"^ for 
a 10"^ K partially ionized medium, to ~ 50 kms"^ for a warm 
10^ K ionized one. Given the typical values of proper motion, 
a pulsar drives a bow-shock in the ISM, that might be revealed 
as a Ha nebula, with emission mainly due to charge-exchange 
and collisional excitation of neutral hydrogen in the bow shock 
itself. 

Only five pulsar wind bow- shock nebulae have been ob- 
served in Ua so far: PSR 1957+20 (Kulkarni & HesterJ988J . 
PSR 2224-h65, the Guita r Nebu la (Cordes et al. fT993l . PSR 
J0437 +4715 (Bell et al. Il995li . PSR 0740-28 (Jones et al. 
I2OOTI1 . and PSR 2124-3358 (Gaensler et al. 2002^. In addi- 
tion, it is not clear yet whether the cometary nebula associated 
with the isolated neutron star RX Jl 85635-3754 (van Kerkwijk 
& Kulkarni 12001ft results from a bow- shock. Finally, the re- 
cently detected small Ha nebula coincident with the position 
of AX J085 1.9-47 17. 4, loca ted at the center of SNR G266.2- 
1.2 (Pellizzoni et al. "2003V has also been interpreted as the 
signature of a new bow- shock nebula. 

In analogy with plerions, one might in principle expect 
bow- shock PWNe to shine also at radio and X-ray frequencies 



2 



N. Bucciantini, E.Amato, L.Del Zanna: RMHD simulations of pulsar bow-shock nebulae 



through synchrotron radiation. X-ray emitting electrons are a 
powerful probe of the fluid dynamics because their short life- 
time allows one to trace the flow structure inside the nebula and 
to constrain the pulsar properties. However Ha bow- shocks are 
usually associated with old pulsars, whose luminosity has of- 
ten fallen below the threshold for which the associated PWN 
could be detected in X-rays. It is thus not surprising that sig- 
nificant advances in detecting non-thermal emission from such 
nebulae have been made after the launch of the Chandra X-ray 
Observatory. 

Some interesting examples of elongated non-thermal nebu- 
lae associated with pulsars still inside the SNR or just emerg- 
ing have been recently found. IC443 (Bocchino & Bykov 2001 ) 
shows an elongated shape and a map of spectral index suggests 
it results from a pulsar moving in its SNR. There has been no 
pulsar detection so far, even if there is evidence for an excess 
of thermal X-ray emission that can be interpreted as the sig- 
nature of a neutron star. Another example of a more distorted 
nebula is W44 (Petre et al. 12002b . where the pulsar is at the 
tip of a narrow protuberance connecting it to the main body 
of the PWN. In CTB80 the pulsar PSR B 195 1+ 32 is located 
just inside the outer edge of the SNR (Fesen et al. 1988V There 
is evidence, from synchrotron emission, for a PWN inside the 
remnant, but the position of the pulsar lays just on the edge of 
a narrow emission finger that protrudes from the main body of 
the PWN. This was interpreted assuming a higher speed and 
the interaction with the denser SNR shell material. 

A possible transition case between pulsars inside and out- 
side their SNR is that of G5.4-1.2, although recent measure- 
ments of the pulsar proper motion have questioned the PSR- 
SNR association. Here the pulsar, that appears to have just 
emerged from the SNR, is at the tip of a flat spectrum radio pro- 
tuberance just outside the SNR shell (Frail & Kulkarni fl991ll . 
Coincident with such radio finger there is a X-ray tail that is 
likely to be synchrotron radiation from the shocked pulsar wind 
(Kaspi et al. 1200 111 . The X-ray nebula does not seem less ex- 
tended than the radio tail, thus suggesting that the fluid speed 
must be high enough so that the flow time is shorter than the 
high energy electrons lifetime. 

The first detection of X-ray emission from a bow-shock 
nebula interacting with the ISM was associated with the Guitar 
Nebula (Romani et al. fT997l based on ROSAT data). However 
recent Chandra observations have shown the presence of a faint 
X-ray filament in the field (Chatterjee 2004) which is clearly 
not associated with the nebula, and which might be at the origin 
of that first detection. Recently, X-ray emission has been de- 
tected in the vicinity of PSRB1957+20 (Stappers et al. 2003), 
a pulsar already known to be associated with a Ha bow-shock 
nebula. There is a bright extended X-ray source coincident with 
the position of the pulsar itself and this has been interpreted 
as the signature of the wind termination shock (TS). A X-ray 
tail aligned to the bow- shock axis is seen to extend from this 
nebula in the direction opposite to the pulsar proper motion. 
However, perhaps the most impressive non thermal bow-shock 
nebula, due to the high quality of both radio and X-ray data, is 
G359.23-0.82, the Mouse Nebula (Gaensler et al. 2004). Again 
the X-ray tail is very extended compared to the assumed size 
of the bow-shock. The X-ray emission peaks close to the pul- 



sar, where the TS is expected to be. Radio imaging of the tail 
seems also to suggest a structured morphology with a fainter 
and larger envelope. 

In all these cases it seems that the X-ray tail has an almost 
cylindrical shape and a smaller transverse section than what 
expected based on analytic models of bow- shocks. Moreover 
observations seem to suggest that the flow velocity in the tail 
must be much higher than the proper motion of the pulsar, and 
possibly close enough to the speed of light to give appreciable 
beaming. 

A difl'erent case is that of Geminga, where two faint tails of 
emission, extending in a direction opposite to the pulsar proper 
motion, have been revealed by Chandra (Caraveo et al. 2004 ). 
This emission might be connected with a possible bow-shock 
nebula, despite the fact that Geminga is not a typical radio pul- 
sar, and that the two tail morphology is difl'erent from what is 
observed in other bow-shock nebulae both outside and inside 
the SNR. 

Apart from PSR B 1957+20, in no other case non-thermal 
and Ha emission have been detected together. This dichotomy 
might be explained considering the efficiency for producing 
these two radiation components and how it evolves with the age 
of the pulsar. Young pulsars have a higher E, favoring the de- 
tection of non-thermal emission, but they have a higher surface 
temperature (Tsuruta"1986), and the photo-ionizing flux is usu- 
ally high enough to completely pre-ionize the ISM in front of 
the bow- shock. Older pulsars are usually much colder, so that 
neutrals can survive (Bucciantini & Bandiera 2001 ), but they 
are characterized by a lower value of E (we are not considering 
here millisecond pulsars that belong to binary systems) so that 
even if a synchrotron nebula is present the emission easily falls 
below the threshold of detectability. 

Although it is known that pulsar winds are magnetized, and 
that the ratio between the Poynting and kinetic energy flux in 
the wind can in principle be as high as 1 (see i.e. Arons 20021 
for a discussion on pulsar wind magnetization and the cr prob- 
lem), models for these nebulae have been presented only in the 
classical HD regime (Bucciantini 2002a, van der Swaluw et 
al. I2003II . This can be considered a good approximation for 
the outer layer of shocked ISM, where the Ha emission comes 
from: here one does not expect the magnetic field to play a 
dominant role given the typical low magnetization of the ISM. 
However the neglect of magnetic fields, as well as the non rela- 
tivistic assumption, in principle, prevent one to use those mod- 
els to infer the properties of the relativistic shocked wind, and 
to obtain more realistic results to compare with observations, 
as far as non-thermal emission is concerned. 

The general structure of the fluid dynamics in pulsar bow- 
shock nebulae is a 3D problem. There are two preferential di- 
rections: the pulsar velocity direction, which defines the "sym- 
metry" axis of the cometary nebula; and the spin axis of 
the pulsar, which defines the direction of the wound-up mag- 
netic field. Additional 3D eff'ects may arise also as a con- 
sequence of anisotropic pulsar wind energy flux (Bogovalov 
fT9^ Komissarov & Lyubarsky'^OO^^ Del Zanna et al. 2004). 
The complete treatment of the system in 3D is at present 
computationally prohibitive. The difl&culty is intrinsically con- 
nected with the presence of two very different time scales. The 
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time scale for the relativistic wind component, which defines 
the time step, is much smaller than that of the non relativistic 
ambient medium, which defines the global evolutionary time 
scale of the nebula. A simplification is possible if one assumes 
that the spin axis is aligned with the pulsar velocity direction. 
In such case the problem is reduced to an axisymmetric geom- 
etry which can be easily modeled in 2D. 

We want to stress that such an approximation, in addition to 
being supported by observations of the Vela and Crab pulsars 
(Helfand et al. 2001), can be also justified on the basis of the 
models that until now have been proposed for the origin of the 
pulsar proper motion. These models can be essentially divided 
in three categories, depending on the origin of the kick: super- 
nova explosion asymmetries (Lai & Goldreich 2000), disrup- 
tion of a binary system (Tauris & Takens[l998), MHD rocket 
eff'ect (Harrison & Tademaru l 1975ft . In general, if the source of 
the kick lasts longer than the core collapse (1 ms) one expects 
the spin axis to end up parallel to the direction of motion of the 
pulsar itself (Spruit & Phinnev fT998t . 

In the following we present the results of a series of 2D 
relativistic MHD simulations of these systems, focusing on the 
eff'ects of the degree of magnetization of the pulsar wind. In 
Sect. 2 we describe the setup and the initial conditions of our 
numerical models, as well as discussing the limitations intrin- 
sic to our approach. In sect. 3 we discuss the changes in the 
flow structure and in the behavior of the various fluid quan- 
tities for increasing values of the magnetization parameter. In 
Sect. 4 we compute simulated emission maps resulting from the 
flow structure we find. We discuss the general spatial proper- 
ties of the synchrotron emission in these systems and how the 
extent of the nebula, is related to the wind magnetization. A 
comparison with observations in the case of the Mouse Nebula 
is presented. Finally, in Sect. 5 we summarize our conclusions. 

2. Simulation setup and initial conditions 

Simulations are performed using the shock-capturing cen- 
tral scheme for relativistic MHD developed by Del Zanna & 
Bucciantini (2002| and Del Zanna et al. (I2003II . to which the 
reader is referred for the numerical algorithms employed. If 
one assumes the pulsar spin axis to be parallel to the pulsar di- 
rection of motion, the flow dynamics can be reduced to a 2D 
axisymmetric configuration in the frame comoving with the 
pulsar. In this case one can use a cylindrical grid (R,z,(p). A 
further simplification of the problem can be achieved if one as- 
sumes V(i) = Bz = Br = 0. In this case the number of fluid 
variables evolved by the code is reduced from 8 to 5, and the 
fluid equations are simplified (see also Del Zanna et al. 2004). 
Symmetric reflection conditions are imposed on the axis, while 
continuous zeroth-order extrapolation is imposed on the other 
boundaries of the computational box. We also assume both for 
the pulsar wind and ambient medium the same adiabatic index 
4/3 (see Bucciantini 2002a and Bucciantini et al. 2004 for mul- 
tifluid simulations with two difl'erent adiabatic coefficients and 
the related problems). 

A spherically symmetric relativistic wind is injected from 
location (R = 0,z = 0) representing the pulsar position, 
with the following characteristics. The wind Lorentz factor is 
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7=10; the wind is cold, with pc^/p = 100. In order to evaluate 
how the flow structure changes with the wind magnetization, 
three values of the magnetization parameter cr = r^B(r)^c/E 
are adopted, namely cr = 0.2, 0.02, 0.002. The toroidal field, 
wrapped around the symmetry axis, can be assumed to have 
the same polarization on both hemispheres since the MHD 
equations, in our simplified 2D axisymmetric setting, only de- 
pend on B^. The ambient medium is modeled as a uniform un- 
magnetized wind entering our computational box with a speed 
V = 0.03c and Mach number M = 30. Even if the wind speed 
is about one order of magnitude higher than the proper motion 
of the fastest known pulsar, the outer medium still evolves as a 
non-relativistic fluid. This choice is only made to allow a faster 
evolution of the nebula without changing significantly the flow 
dynamics with respect to a case where a typical value of the 
pulsar speed is used. 

As it is known both from analytic theory (Wilkin 119961) 
and from classical HD simulations rBucciantini l2002ai van der 
Swaluw 2003 ), bow-shock nebulae are characterized by a typi- 
cal length scale, the so called stagnation point distance, dg. This 
is the distance from the pulsar at which the wind momentum 
flux is balanced by the ISM ram pressure: 

do = ^EliAncpoV^) (1) 

with po the density of the ISM. This means that the dimensions 
of the nebula are a function of the total energy flux from the 
pulsar, and the ram pressure of the ambient medium only. We 
have chosen a grid with uniform resolution (AT? = Az) and a 
number of cells such that J^/Az = 80. In a HD case this cor- 
responds to having about 100 cells along the axis between the 
pulsar and the contact discontinuity (CD) position (Bucciantini 
I2002al . whereas this number of cells increases with cr as we 
will discuss in the next section. 

The axisymmetric assumption, which provides an essential 
simplification of the problem, forces the magnetic field to be 
exclusively toroidal and wrapped around the symmetry axis. 
The reduction of the degrees of freedom of the system leads 
to unphysical over- stability, that has important implications on 
the dynamics especially in the head of the nebula. Here, due 
to hoop stresses, magnetic field loops tend to pile up. For low 
magnetizations this has the eff'ect of unrealistically deforming 
the very head of the system, close to the axis, into a cone-like 
shape. With increasing values of cr, no "quasi- steady" solution 
can be found, and the CD is pushed out of the computational 
box. We considered increasingly large computational domains, 
up to the point of placing the left edge of the computational 
box as far as 25do from the pulsar in the case a = 0.02. This 
fact did not solve the problem but allowed us to measure the 
speed at which the CD moves out of the box at several diff'erent 
times. We find that this speed (which is a fraction of the speed 
of light) does not decrease in time, contrary to what one would 
expect if a steady state were to be reached at some point. We 
interpret this fact as a result of magnetic collimation: when the 
latter becomes important the post-shock flow in the vicinities 
of the axis is not allowed to diverge laterally contrary to the 
requirements for the existence of a steady state solution. 
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In a real multidimensional case, kink instabilities will even- 
tually arise (Begelman 1998), destroying the symmetry of the 
problem, partly dissipating magnetic field. In our simplified ge- 
ometry, where this kind of instability is not allowed to grow, 
we decided to artificially dissipate magnetic fields. Our choice 
is justified by the expectation that kink instability will eventu- 
ally arise in a strong magnetic filed, however the dissipation 
we introduced has no pretence of reproducing the realistic de- 
velopment of kink instability. In the case with cr = 0.002 it is 
possible to find a quasi-steady state solution without imposing 
any dissipation. Whereas for the case with cr = 0.2 we had 
to modify the wind structure by introducing an unmagnetized 
conical region: within an angle of 20 deg from the axis in the 
forward direction the wind is unmagnetized, though its energy 
flux is the same as everywhere else in the solid angle; moreover 
the magnetic field strength is set to zero also downstream of the 
TS, mimicking magnetic dissipation. To allow direct compari- 
son among the various cases we have investigated, we decided 
to adopt the modified wind structure for all values of cr. 

3. Flow structure 

The geometry and internal flow structure of the bow- shock neb- 
ula are shown in Figs. lllSl for the three values of the wind mag- 
netization we adopted, namely cr = 0.002, 0.02, 0.2. 

3.1. Outer layer 

First of all let us discuss the shape of the external layer, which 
is where Ha emission is believed to originate. We notice that 
this does not change significantly with the wind magnetiza- 
tion. Major deviations are present only in the head where hoop 
stresses and magnetic pinching are stronger. Without magnetic 
dissipation, the piling up of magnetic loops drives the front of 
the bow- shock to move further away from the pulsar along the 
axis, causing the very head to assume a rather conical shape. 
When artificial dissipation is included, the global shape (es- 
pecially in the tail) does not show sensible variations as the 
wind magnetization increases. However the front shock dis- 
tance from the pulsar (along the axis) can change significantly. 

In the HD case, it can be proved that, if both the ambient 
medium and the stellar wind have the same adiabatic coeffi- 
cient, the stagnation point distance Eq. [I] exactly defines the 
distance of the wind TS from the pulsar along the axis in the 
forward direction (FTS). However, even if the adiabatic coef- 
ficients of the two media are diff'erent, JpTS - still holds 
as an approximate relation. The positions of the CD and of 
the front bow- shock depend on the fluid properties (i.e. com- 
pressibility) of the wind and ambient medium respectively. An 
ordered magnetic field is less compressible than a gas and its 
pressure, rather than increasing above equipartition, pushes the 
front bow- shock further away from the pulsar as cr increases. 

It is interesting to notice that there is a "critical" polar angle 
for hoop stresses: streamlines originating closer to the axis are 
collimated toward it while those originating at larger polar an- 
gles are advected in the tail. We found that such critical angle 
increases with magnetization and its value is 5°, 20° and 80° 
for cr = 0.002, 0.02, 0.2 respectively. 



The increase of the distance of the front shock from the pul- 
sar has consequences that go beyond changing the length- scale 
of the system. It aff'ects both the microphysical properties of the 
outer layer and its Ha emission, and the inferred conditions of 
the ambient medium. While the pulsar spin down energy and 
its proper motion can be measured directly, the front shock dis- 
tance is used to infer the ISM density. But increasing the mag- 
netization of the wind from (T = 0tO(T = 0.2, we find that such 
distance increases by a factor ~ 3: consequently the inferred 
outer density might be underestimated by about one order of 
magnitude, if it is computed based on the HD model. The same 
kind of mistake could aff'ect the estimate of the density of the 
neutral component based on Ha emission, given that the sur- 
face density of the external layer, which determines the fraction 
of interacting neutrals, scales as podo rBucciantini l2002bll . 

3.2. Pulsar wind cavity 

Let us focus now on the internal portion of the nebula occu- 
pied by the magnetized relativistically hot plasma coming from 
the pulsar wind. It has been suggested that the shape of the 
TS can be observed in the X-rays. For instance, Stappers et al. 
l2003l interDret the luminous knot around PSRB 1957+20 as the 
TS itself, while in the case of the Mouse Nebula Gaensler et 
al. (I2004II suggest the identification of the TS with the X-ray 
tongue. This latter work shows that the elongation of the X- 
ray emitting region agrees with the shape of the TS one could 
expect based on a classical HD model. In previous numerical 
studies, Bucciantini (2002a) and van der Swaluw et al. (2003]) 
assumed that asymptotically the pressure in the tail would reach 
the value of the ISM pressure, and found the elongation of the 
TS (ratio between the distance of the forward and backward 
parts of the TS) to be proportional to the Mach number M. 
In the simulation by Gaensler et al. (2004), instead, the back- 
ward TS (BTS) never moves further from the pulsar than 5-6 
times the FTS distance, for whatever value of the Mach num- 
ber. Our simulations confirm that the average pressure in the 
tail seems to saturate at a value Ptaii ~ 0.02poV^. The detailed 
pressure profiles in the z-direction show a shallow gradient: this 
leaves open the possibility that pressure equilibrium with the 
ISM is reached, at large distances. However the elongation of 
the TS seems to depend strongly on the magnetization of the 
nebula. With increasing magnetization the BTS moves toward 
the pulsar and the back of the free flowing wind cavity shrinks. 
Concerning then the FTS, in our simulations this keeps approx- 
imately the same distance from the pulsar for all values of cr. 
However the efl'ect of the artificial dissipation we introduced 
is certainly important for the fluid dynamics in the head, hence 
no definite conclusion can be drawn on the position of the FTS, 
although we expect it to be less aff'ected by the magnetization, 
than the BTS, since, as stressed before, its position is deter- 
mined mainly by ram pressure equilibrium. 

Also the shape of the BTS changes with magnetization. 
In HD cases it is concave toward the pulsar and the shock is 
perpendicular to the wind speed (Bucciantini 2002a). As cr in- 
creases the shock becomes convex. This is a consequence of 
the radial profile of the total pressure in the nebula. When the 
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azimuthal magnetic field is present the internal pressure is non- 
uniform in the R direction, tending to be higher toward the axis 
(e.g. Begelman & Li 1 1992ft : this causes the BTS on the axis to 
be closer to the pulsar than at larger cylindrical radii. As the 
BTS becomes more and more convex it also becomes highly 
oblique and the parallel (to the shock) component of the veloc- 
ity becomes a non negligible fraction of c. 

3.3. Inner flow structure and tail 

Regarding the general shape of the inner region occupied by 
the hot relativistic plasma, these new simulations confirm pre- 
vious results (Bucciantini 2002a) suggesting that the material 
is collimated in a cylindrical tail. The radial extent of this cylin- 
der is found to be 7?taii - 4(io and it does not depend neither on 
the wind magnetization nor on the Mach number (in the hyper- 
sonic limit). This can be easily understood in terms of energy 
conservation and general scalings for the fluid quantities. In the 
weakly magnetized cases one can write the energy conservation 
as: 



E = 27T 



f 

Jo 



(2) 



with all quantities in the integral depending in general on the 
cylindrical radius R, Vtaii being the flow speed. Neglecting the 
radial profiles and using the definition of the stagnation point 
(Eq.ID: 



PoV^ 



1 c 



dp 



(3) 



The only scale velocity for the plasma inside 7?taii is the speed 
of the wind which is ~ c. The "two thin layers" solution 
(Bucciantini I2002bl would predict Vtaii ~ 0.6c. In fact, in our 
simulations (middle panel of Figs. 1113b we have, in all cases, 
Vtaii ~ 0.8 - 0.9c (ytaii ~ 2). In addition we find, as mentioned 
above, Ptaii ~ 0.02poV^. Substituting these numerical values in 
Eq.lSjone easily obtains 7?taii ~ "^dg. 

Considering in greater detail the internal structure of the 
relativistic plasma flow in the tail, we find a qualitative agree- 
ment with previous analyses (Bucciantini 2002a). The flow has 
difl'erent average properties in two main regions (see middle 
panel of Figs. lllSll . The inner channel contains particles com- 
ing from the BTS. The latter, in the HD and cr = 0.002 cases, is 
almost perpendicular to the wind propagation direction, so that 
the flow downstream is subsonic (flow speed ~ c/3 in the HD 
case). The outer channel, surrounding the former, contains par- 
ticles coming from the head of the nebula and from a portion 
of the TS which is oblique with respect to the wind speed. The 
fluid in the head is subsonic due to the shock geometry, but, as 
it expands sideways, it may undergo a sonic transition and be- 
come supersonic, or superfastmagneto sonic in the MHD case 
(sonic surfaces are represented as contours in the upper panel 
of Figs. fTlSl . The angular distance from the apex at which this 
happens depends on the magnetization. For low cr it takes place 
at a polar angle ~ 60° - 70°. As the magnetization increases 
the transition occurs further away: for cr = 0.02 the angle is 
~ 120°, while for cr = 0.2 the downstream flow is everywhere 
subfastmagnetosonic and transitions take place in the tail. 



Typical speeds of the two channels are difl'erent. In the su- 
personic outer channel the velocity is not much dependent on 
magnetization and is of order 0.8 - 0.9c. In the subsonic in- 
ner channel values of the speed are in the range 0.1 - 0.3c 
for cr = 0.002 and in the range 0.2 - 0.5c for cr = 0.2. The 
shear between the flow in these two channels is at the origin of 
Kelvin-Helmoltz type instabilities at their CD. The efl'ects of 
this instability are not evident in the top panels of Figs. fTl3l due 
to the low density contrast, but with a greater contrast the pres- 
ence of mixing and dragging from the fast component to the 
slow one becomes clear. Moving away from the TS, the speed 
in the slow channel tends to increase as a consequence of this 
dragging and it may become supersonic. 

Shear instability is also present at the CD between the outer 
layer of shocked ISM and the inner layer of shocked pulsar 
wind. As we mentioned, the flow speed in the latter is close to 
c, while in the former the fluid velocity is of the same order 
as the pulsar speed through the ISM, V. The large velocity dif- 
ference makes the development of Kelvin-Helmoltz instability 
at the interface an obvious expectation. However no signs of it 
had been observed in previous numerical studies (Bucciantini 
|2p02a van der Swaluw 2003), where the flow in the tail was 
completely laminar. This may partly be due to the large artifi- 
cial and numerical diff'usion present in those works, and partly 
be intrinsic to the non-relativistic treatment of the problem. In 
fact, in the framework of classical HD, the velocity diff'erence 
at the CD is smaller and the density jump between the two lay- 
ers is strongly reduced. 

In a proper relativistic treatment the ratio between the en- 
thalpy of the outer (~ p^c^) and inner layer (~ 4P ^ 4poV^) is 
~ v^/c^. For typical values of pulsar proper velocities V ~ 300 
kms"khis ratio is 10^. This implies that even a very weak 
contamination (1 part in a million) of the inner region by the 
shocked ISM might substantially aff'ect the internal flow dy- 
namics, while larger contamination can in principle modify the 
overall nebular shape (Lyutikov 2003 ). In our case fingers and 
eddies of heavier material (that in the simulation is however 
about 1000 times less dense than realistic ISM) can detach from 
the outer layer forming clumps of non relativistic material em- 
bedded in the hot magnetized pulsar wind plasma. Simulations 
show that the instability forms at intermediate latitudes in the 
portion of the head of the nebula that is still subsonic. In the 
very head, close to the apex, the shear velocity is smaller and 
the instability is not present, while in the supersonic part of 
the tail clumps of material are advected far from the pulsar at 
a speed ~ 0.5c. These clumps can in principle aff'ect the flow 
structure in the tail while they are advected away: this effect 
could show as a modulation in the synchrotron emission maps. 

Shear instability is also at the base of two other phenom- 
ena observed in our simulations. First of all, it is responsible 
for the fluctuations in the shape of the TS. Secondly, it is at the 
origin of a series of oblique shocks that propagate in the outer 
supersonic channel. Such oblique shocks are a common feature 
in supersonic collimated flows, and are found in many jet sim- 
ulations (see for example Komissarov 1999). Even if they are 
weak or only moderately strong, they are important as possible 
sites of particle acceleration. 
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Fig. 1. Flow structure in the case cr = 0.002. Upper panel: the 
color map represents the density in logarithmic scale (arbitrary 
units); the contours define the sonic surface. Middle panel: the 
color scale refers to the velocity magnitude; the arrows trace 
the streamlines. Bottom panel: map of the magnetization in the 
nebula defined as B^/i^ny^P). Spatial scale is in units of 1.3 do, 
corresponding to the distance of the CD from the pulsar, along 
the axis, in the HD case. 



However, we must point out that the setting of our simu- 
lations does not guarantee a proper treatment of either of the 
shear unstable surfaces we see, both between fast and slow 
channel and between the outer and inner layer. The main limi- 
tation arises from the 2D axisymmetric approximation and the 
absence of a poloidal component of magnetic field that could 
stabilize the interfaces. We also observe that the instability in- 
creases with magnetization, this might in principle be due to the 
fact that for high cr the flow is subsonic in a more extended re- 



Fig. 2. Flow structure in the case cr = 0.02. The legend is the 
same as for Fig. [l] Now in the bottom panel contours appear, 
that represent equipartition surfaces. 



gion, although, once again, numerical problems cannot be com- 
pletely excluded as a cause. 

3.4. Magnetic field distribution 

Finally, let us consider the structure of the magnetic field. As it 
is known from solar wind simulations, the magnetic field at the 
CD with the outer layer is compressed. An interesting ques- 
tion is how thick this magnetic sheet is. If the magnetization 
in the wind is not too high (cr < 0.1) the flow structure is es- 
sentially the same as in a HD case. Magnetic field is passively 
advected and compressed toward the CD where it eventually 
reaches equipartition. As one can see from the bottom panel of 
Fig. [l] equipartition is reached in a narrow sheet close to the 
CD. In the lowly magnetized cases the nebular magnetic field 
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is on average well below equipartition. This behavior changes 
when cr is larger than ~ 0.1. This corresponds to the value at 
which the magnetic field in the tail can be estimated to be at 
equipartition. From flux conservation the magnetic field in the 
tail of the nebula, Btaii, is, in the pulsar frame: 
r^w(r) c 

^tail = 71— (4) 

^tail Vtail 

with Bw(r) the magnetic field in the wind at a distance r from 
the pulsar {B^{r) oc r"^), and Rtaii the radial extent of the tail, 
^taii - ^do. When, in the tail, the thermal pressure, Ptaii, is 
dominant, its average value can be derived from energy conser- 
vation, Eq.|3l Neglecting once again the radial profiles: 



Density (LoglO); a = 0.2 



P tail - 



We then obtain for the magnetization in the tail: 



^tlii 



6cr. 



(5) 



(6) 



;ail 



When cr ^ 0.1, the dynamics in the tail is dominated by the 
magnetic field. Pressure in the internal region is no longer 
nearly uniform in the transverse direction, and becomes larger 
toward the axis, due to magnetic pinching, causing the BTS 
to recede and become convex in shape. Moreover, while for 
cr = 0.02, 0.002, equipartition is reached only in the vicinity of 
the CD, now the magnetic field is above the equipartition value 
also close to the axis. Since synchrotron radiation depends on 
B^, this change in the magnetic energy profile might give im- 
portant signatures in the emission maps. 

4, Synchrotron emissivity 

Given the flow structure, magnetic energy distribution and 
pressure profile inside the nebula, it is possible to compute syn- 
chrotron emission maps. Following Begelman & Li (1992), we 
assume that particle acceleration takes place only at the TS and 
the resulting particle energy distribution is: 

with 

= 3Kp(2a-l)- 



(7) 



for a it 0,5 . 



(8) 



Here Kp PoC^) is the fraction of the thermal pressure down- 
stream of the TS that is carried by accelerated particles, 6^ and 
6m are the minimum and maximum energy of the particle spec- 
trum respectively, and ^ labels the diff'erent streamlines diverg- 
ing from the TS. Po is diff'erent all along the TS but since the 
shock is everywhere strong, we assume^^, and em to be in- 
dependent on the streamline. Moreover in the following we will 
always assume <^ €m and a 0.5. 

At a generic distance s along a streamline the particle dis- 
tribution function, evolved taking into account adiabatic and 
synchrotron losses, is : 




12 3 4 

Velocity magnitude + streamlines; a = 0.2 
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Fig. 3. Flow structure in the case cr = 0.2. The legend is the 
same as for Figs. 11121 



with 



6oo = 6oo(^,^)=p'/^ 



4 e 

9 m'^c 



ds 



(10) 



where p is the relativistic plasma density. As we will discuss be- 
low, for the typical range of pulsar velocities and luminosities, 
synchrotron losses are negligible within the spatial region cor- 
responding to our computational domain and for the frequency 
range we are interested in. This fact translates into Eoo » 6, that 
allows one to simplify Eq.|51 

Another eff'ect that needs to be taken into account is 
Doppler boosting, since the emitting plasma is moving with 
bulk velocity close to c. We decided for simplicity to consider 
a nebula being viewed face on, i.e. with the pulsar speed lay- 
ing in the plane of the sky. In calculating the volume emissivity 
s{v) toward the observer we have approximated the radiation 
spectrum of the single particles as a Dirac delta centered on the 
critical frequency Vc = C2Be^. In this approximation we obtain: 



X (P(^,^)/W))—(l- 6/600) 



2a-l 



(9) s{y) 



ci c: 



a-l 



a+l 



(11) 
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where ci = (4/9)e'^c/(mc^f, C2 = 0.29 x 3ec/(47i(mc^)^), B'^ is 
the magnetic field component perpendicular to the observer's 
line of sight, as measured in the fluid frame, and Z) = (7(1 - n • 
is the Doppler factor. Here P is the flow velocity and n is 
the unit vector pointing toward the observer, both in the pulsar 
frame. The latter is related to the observer direction in the fluid 
frame through (see e.g. Komissarov & Lvubarskv l2QQ4ll : 



R2 



i)-r 



(12) 



Before computing the emission map at a given frequency 
we need to choose a spectral index a. Estimates of the X-ray 
spectral index are presently available only for four pulsar bow 
shock nebulae: for IC443, Bocchino & Bykov (2003 ) find a < 
0.5; for the Duck nebula Kaspi et al. dJoHD give a = 0. ± 0.6; 
for PSR B1957+20 Stappers et al. (tM^ find a = Q.9 ± 0.5; 
finally, Gansler et al. (2004) give for the Mouse nebula (0.8 ± 
0.1) <a< (1.6 ±0.2). For our simulated map we decided to use 
a = 0.6, a choice which avoids the critical value of = 0.5 and 
leaves in Eq.[n]only a very weak dependence on the conditions 
at the TS (PJ. 

The synchrotron emission maps, that we built by integrat- 
ing the volume emissivity in Eq.[n] along the line of sight, are 
shown in Fig.|4]for all the diff'erent values of cr we adopted. The 
brightest part of the nebula corresponds to the head of the bow- 
shock, where magnetic field and particle pressure are highest. 
Here Doppler boosting is not very important since the flow ve- 
locities involved are low. Moving from the head toward the 
back of the nebula we find a region of fainter emission and then 
an even fainter tail. It is interesting to notice that in the case 
of the Mouse nebula Gaensler et al. (2004) distinguish three 
main regions in the X-ray emission pattern, that they name 
as the "head", the "tongue" and the "tail", with intensity ra- 
tio 100/10/1. This is the same intensity ratio we find in our 
synthetic maps between the corresponding regions. The ratio 
between the head and the tongue, moreover, turns out to be the 
same for all values of cr. This would imply that the head/tongue 
emission properties do not provide us with information on the 
wind magnetization. 

We also notice that the tail of emission in the back of the 
nebula is more enhanced toward the axis than close to the CD. 
This is a consequence of the toroidal structure of the field. 
Close to the axis the magnetic field lines are perpendicular 
to the observer and (5^) is maximum, while at the edge of 
the inner region the magnetic field is almost pointing toward 
the observer, so that (5^) reaches a minimum. However, as 
we stressed in the previous section, shear instabilities both at 
the CD and between the fast and slow internal channels can in 
principle destroy the ordered structure of the magnetic field. If 
a disordered component of the magnetic field is introduced in 
Eq. [n] the contrast between the axis and the edges weakens. 
In Fig. 13 we show the synchrotron map corresponding to the 
case cr = 0.002 when a completely disordered magnetic field is 
used in Eq.[n] while keeping the magnetic energy distribution 
inside the nebula the same. 

The diff'erence in the axis-edge intensity contrast between 
an ordered and a disordered magnetic field becomes smaller 
with increasing cr. This is because, as magnetic stresses in the 



tail become dynamically more important, the magnetic energy 
density becomes higher close to the axis, where projection ef- 
fects are less important. For the same reason, in comparison 
with cases at lower cr, for cr = 0.2 the emission is enhanced 
toward the axis (bottom panel of Fig. |4|). We suggest that the 
observation of an emission tail narrower than the head of the 
nebula might be considered a signature of high (> 0. 1) cr winds. 
However a high cr wind is not bounded to produce such a struc- 
ture, because this is strictly related to the axisymmetric approx- 
imation. In fact, if the spin axis is not aligned with the pulsar 
velocity, the magnetic energy distribution is going to be dif- 
ferent. Moreover, even in the case of alignment, it is always 
a possibility that 3D instabilities disrupt the ordered magnetic 
field structure eliminating hoop stresses. 

As we discussed in the previous section, the shear at the CD 
is at the origin of instabilities and causes the formation of dense 
clumps of unmagnetized shocked ISM, that are embedded in 
the relativistically hot plasma. Such instabilities, that appear to 
be stronger at high magnetization, and the clumps they gener- 
ate, can in principle modulate the emission in the tail, both dis- 
torting the fluid flow pattern and compressing the relativistic 
material. Typical time scales for the detachment of clumps and 
their flow time in the tail are of order dole. Evidence of modu- 
lation in the emission pattern of non thermal radiation in pulsar 
bow- shock nebulae has been found in the case of the Mouse. 
Unfortunately we have no information about the magnetic field 
structure in this nebula (i.e. pulsar spin axis inclination) and/or 
time variability of the observed modulation that could confirm 
or rule out theories on its origin. In fact other possible explana- 
tions have been suggested so far that do not require instability 
at the CD, as expansion-compression waves (Lyutikov 2004 ). 

As a final remark on the spatial structure of the emission let 
us shortly discuss the dependence of the results on our assump- 
tion for the particle spectrum and for the orientation of the neb- 
ula. In the case of a particle distribution function with a much 
diff'erent from 0.5, in building emission map the conditions at 
the TS would play a much more relevant role (see Eq.[n}. In 
general the eff'ect is to enhance the contrast between the vicini- 
ties of the pulsar and the edges of the nebula. However, even in 
the case of particles coming from the head, where the pressure 
is highest, and flowing in the tail, the correction is just a factor 
~ 2 for Of = 1 . Much more dramatic changes in the emission 
pattern may come instead from a difl'erent orientation of the 
pulsar velocity. If the pulsar proper motion has a non-negligible 
component out of the plane of the sky, relative motions in the 
tail become important. For a pulsar moving toward the observer 
the emission in the tail is weaker than what is shown in Fig.|4l 
while in the case of a pulsar moving away from the observer 
emission should be enhanced. 

Let us focus now on the issue of synchrotron losses. In the 
above discussion we have assumed that particles of all ener- 
gies, within the range of interest for emission in the X-rays, 
survive everywhere in the computational domain. An estimate 
of the synchrotron lifetime might be easily given. The magnetic 
pressure in the the tail of the nebula is approximately given by 
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Eq.|6| and, using Eq.[3l can be written as: 



i^Vtaii; \^tail/ 



\Vtail/ 



(13) 



We approximate the age of a particle found at a distance rj^do 
from the apex as ~ i^d^o) / (^v<^) where rj^c is the average flow 
speed the particle has experienced. In the following we will 
be interested only in particles flowing in the tail and with and 
average velocity such that rjy ^ Vtaii/c: this implies that we are 
not considering particles accelerated at the FTS, within a polar 
angle of about n/4 from the axis. Given the magnetic field in 
Eq. the maximum energy of the photons emitted by this 
particle turns out to be: 



3xl0V-i,-(^) 
\10^ms-W \1034erg/,s-i 



Po 



10-24gcm-3 
1 

keV. 



(14) 



As already pointed out, the value Vtaii - 0.8 - 0.9c is in- 
dependent of the magnetization. Therefore, knowing the extent 
of the tail of the nebula at a given frequency, the above relation 
fEq.nH allows one to estimate the wind magnetization cr. 

4.1. The mouse 

As a special case let us consider the bow- shock nebula around 
PSR J1747-2958 (Gaensler et al. l2004l . also called the Mouse 
Nebula. As we mentioned in the introduction, this is the only 
bow- shock nebula for which data are available both at radio 
and X-ray frequencies. We now wish to discuss these data in 
the light of our numerical results. 

It has been suggested that, in this nebula, synchrotron losses 
are responsible for the diff'erent morphologies observed in these 
two bands. 

First of all, let us consider the most recent X-ray results 
concerning the extent of the nebula. Gaensler et al f2004) give 
for the length of the tail Itaii - 1 . 1 pc, and for the stagnation 
point distance do - 0.024 pc, for an assumed distance of 5 
kpc. Using Eq. [HI with 6^ = 2.5 keV, E = 2.5 x lO^^ergs'^ 
Vo = 6x lO^cms"^ andpo = 0.3 x 10"^'^gcm"^ we find: 



-20(^) 



10/3 



(15) 



As we mentioned above, our simulations give for the flow 
speed in the tail Vtaii/c - 0.8 - 0.9. This implies that the above 
relation cannot be satisfied for cr < 1. It must be stressed, 
though, that Eg.fTSlonlv holds in the low cr limit, /. ^. cr ^ 0.1 
(see Eq. |6|l, and, even using the upper limit, we find that the 
velocity should be Vtaii < 0.2c, in order to obtain the observed 
value for the X-ray extent of the nebula. 

Such a value of the velocity is much lower than what we 
find in our ideal MHD models. A means of slowing down the 
flow efl&ciently is possibly through contamination by the ex- 
ternal non-relativistic ambient medium. This could in principle 



provide at the same time an explanation for another discrep- 
ancy we find between our ideal MHD picture and the data, i.e. 
the radial extent of the tail at radio frequencies. The latter in 
fact is about a factor 3-4 larger than what we find, although pre- 
liminary results of an investigation currently underway show 
that this discrepancy could be alleviated by introduction of a 
latitude dependence of the wind energy flux (analogous to what 
is done for Crab-like plerions, e.g. Del Zanna et al., 2004). The 
flow divergence can also justify the progressive steepening of 
the X-ray spectrum moving away from the pulsar along the 
nebula, as well as the smaller radial extent of the X-ray tail than 
at radio wavelengths. However a self-consistent analysis of the 
efl'ects of mass loading is beyond the scope of this article. 

Using Eqs. lTllll it is also possible to evaluate the expected 
value of Lx/E for the difl'erent values of cr we employed, and 
to compare our results with the value measured for the Mouse 
Nebula, Lx/E = 0.02. We built synthetic maps using the same 
algorithm described above but with a = 1, according to the 
photon index measured by Gaensler et al. ( 2004), and extrap- 
olating the emissivity in the tail to a distance ~ 45 from the 
pulsar. Integration was performed in the Chandra band, from 
0.5 to 5 keV. This allows one to fix the normalization Kp Emin in 
the distribution function and the results are: 



^min 
f^p ^min 



^ W {LxlE)Qxg for a = 0.002, 
^ 10^ (Lx/E)erg for cr = 0.02, 
^5x10^ (Lx/E)QYg for cr = 0.2. 



(16) 



To give an estimate of the fraction of the total energy output 
from the pulsar that must be converted into X-ray emitting 
particles in order to give the observed Lx/E, one can assume 
Emin - 2200 erg, a value corresponding to a particle emitting a 
photon of 1 keV in a magnetic field = ^npoV^. We find that 
even at high magnetization a substantial fraction (~ 30%) of 
the total energy output from the pulsar must be converted into 
X-ray emitting particles. This also implies that the transition 
between the flat (a ~ 0) and steep (a - I) part of the accel- 
erated particle distribution might be close to the X-ray band. 
Again, slowing down of the flow with respect to our findings 
and the consequent compression of the magnetic field might 
make such constraint somewhat weaker. 

5. Conclusion 

In this paper we have presented the first eff'ort to investigate the 
structure of pulsar bow- shock nebulae by means of relativistic 
MHD simulations. Our analysis focused on the signatures of 
the pulsar wind magnetization, cr, on the flow structure and on 
the non-thermal emission by the system. 

We found that the structure of the external layer of shocked 
ISM does not change significantly with varying cr, apart from 
the very head of the nebula. Here magnetic compression is 
stronger and the pile-up of magnetic loops tends to deform the 
rounded HD shape into a pointed one. We can thus conclude 
that the Ha emission should not be much aff'ected by the mag- 
netization of the pulsar wind. While what indeed changes with 
cr is the distance of the front bow- shock, which increases with 
magnetization. This has to be taken into account when using 
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Emission map (LoglO) cj = 0.002 Emission mop (Log10) cj = 0.002 




Fig. 4. Simulated synchrotron maps corresponding to the three 
simulations presented in Sect. 3. The color levels represent the 
intensity in logarithmic scale (arbitrary units). The spatial scale 
is the same as in Figs. I1I3I No disordered component of the 
magnetic field is included in computing the emission. From top 
to bottom: cr = 0.002, 0.02, 0.2. 



pressure equilibrium to evaluate the ambient medium density. 
However, we want to point out that all the effects related to 
magnetic pile-up, especially close to the axis, might be arti- 
ficially amplified in our simulations by the assumption of ax- 
isymmetry, that prevents the development of kink instabilities. 

We found, as expected, that the shape of the TS is sensitive 
to cr. At low cr it tends to be elongated as in the HD case, with 
the distance of the BTS from the pulsar ~ 5 times that of the 
FTS. With increasing magnetization this ratio decreases as a 
consequence of the internal pressure profile. 

In the tail, our simulations confirm the HD result that the 
pressure at the CD saturates at a value ~ 0.02poV^, and the 
expectation that the magnetization increases oc cr. In the low 
cr cases the magnetic field is always below equipartition and 
it tends to peak toward the CD. On the contrary, for cr > 0.1, 
magnetic stresses in the tail are important, the magnetic field 
peaks close to the axis and exceeds equipartition in most of the 
nebular volume. 




Fig. 5. Same as upper panel of Fig. |4] but now the amount of 
magnetic energy that is found at each place in the simulation is 
assumed to be in the form of a completely disordered magnetic 
field. 



The internal flow in the tail is structured in two main re- 
gions: an internal subsonic channel with flow speed ~ 0.2-0.4c 
and an outer fast supersonic channel where velocity reaches 
0.9c. The average flow speed is a substantial fraction of the 
wind speed: we find Vtaii - 0.8 - 0.9c, to be compared with 
the value of 0.6c predicted by the "two thin layers" solution 
rBucciantini l2002bll . 

Emission maps are computed for each flow structure we 
found in our simulations. The main difl'erence when varying cr 
is that, when the wind magnetization is such that most of the 
nebula is above equipartition, the emission is enhanced close 
to the axis. Another point that is worth mentioning is that we 
see modulations of the emission pattern associated with shear 
instabilities in the nebula. In general, the emission maps show 
three main regions, that we denote as a "head", a "tongue" and 
a "tail", in analogy with the labels used for the Mouse nebula 
Nebula. Between these regions the intensity scales, in the order, 
as 100/10/1, a result that is also in qualitative agreement with 
observations of the Mouse Nebula. 

Moreover, the nebular extent at X-ray energies can be used 
to constrain the properties of the plasma in the tail, in terms of 
average speed and magnetization. Comparing our results with 
what is observed in the case of the Mouse Nebula we conclude 
that a low cr wind requires flow velocities in the tail Vtaii ^ c. It 
is not clear how such values can be obtained within the frame- 
work of an ideal MHD model. On the other hand, even as- 
suming a cr that guarantees equipartition in most of the nebula, 
some non ideal process is still required to lower the flow speed 
with respect to the values we find. 

Before concluding, let us stress once more the intrinsic lim- 
itations of our modeling. First of all, the axisymmetric approx- 
imation allowed us to study only configurations in which the 
speed and spin axis of the pulsar are aligned. In addition, the re- 
moval of kink instabilities from the system induces a fictitious 
pile-up of magnetic field that forced us to include "magnetic 
dissipation" in the head in the high cr cases. Finally, we would 
like to remind that for the correct modeling of these systems, a 
non ideal process like mass loading is likely to play a key role. 
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